%matplotlib inline
from scipy import linalg
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML, displaySciPy
Adapted from the notebooks by
Introduction
SciPy is a scientific library that builds upon NumPy. Among others, SciPy deals with:
- Integration (scipy.integrate)
- Optimization (scipy.optimize)
- Interpolation (scipy.interpolate)
- Fourier Transform (scipy.fftpack)
- Signal Processing (scipy.signal)
- Linear Algebra (scipy.linalg)
- Sparse matrices (scipy.sparse)
- Statistics (scipy.stats)
- Image processing (scipy.ndimage)
- IO (input/output) (scipy.io)
Animation display with python
Animation with matplotlib
You can use FuncAnimation to animate a sequence of images:
%%capture
fig, ax = plt.subplots()
xdata, ydata = [], []
(ln,) = plt.plot([], [], "ro")
def init():
ax.set_xlim(0, 2 * np.pi)
ax.set_ylim(-1, 1)
return (ln,)
def update(frame):
xdata.append(frame)
ydata.append(np.sin(frame))
ln.set_data(xdata, ydata)
return (ln,)
ani = FuncAnimation(
fig,
update,
interval=50,
frames=np.linspace(0, 2 * np.pi, 100),
init_func=init,
blit=True,
)display(HTML(ani.to_jshtml()))Another way of displaying video exists, using html5 video:
display(HTML(ani.to_html5_video()))References:
- Matplotlib Animations / JavaScript Widgets by Louis Tiao
Animation with plotly
matplotlib works fine for advanced tuning, but is harder for simple tasks. So just try plotly for basic animations:
import plotly.express as px
from plotly.offline import plot
df = px.data.gapminder()
fig = px.scatter(
df,
x="gdpPercap",
y="lifeExp",
animation_frame="year",
animation_group="country",
size="pop",
color="continent",
hover_name="country",
log_x=True,
size_max=55,
range_x=[100, 100000],
range_y=[25, 90],
)
fig.show("notebook")Linear algebra
scipy for linear algebra : use linalg. It includes functions for solving linear systems, eigenvalues decomposition, SVD, Gaussian elimination (LU, Cholesky), etc.
References:
Solving linear systems:
Find x such that: A x = b for specified matrix A and vector b.
A = np.array([[1, 0, 3], [4, 5, 12], [7, 8, 9]], dtype=float)
b = np.array([[1, 2, 3]], dtype=np.float64).T
print(A, b)
x = linalg.solve(A, b)
print(x, x.shape, b.shape)[[ 1. 0. 3.]
[ 4. 5. 12.]
[ 7. 8. 9.]] [[1.]
[2.]
[3.]]
[[ 0.8 ]
[-0.4 ]
[ 0.06666667]] (3, 1) (3, 1)
Check the result at a given precision (different from ==)
np.allclose(A @ x, b, atol=1e-14, rtol=1e-15)True
Remark: NEVER (or you should really know why) invert a matrix. ALWAYS solve linear systems instead!
Eigenvalues/ Eigenvectors
A v_n = \lambda_n v_n with v_n the n-th eigen vector and \lambda_n the n-th eigen value. The associated python functions are eigvals and eig:
A = np.random.randn(3, 3)
A = A + A.T
evals, evecs = linalg.eig(A)
print(evals, "\n ------\n", evecs)
np.allclose(A, evecs @ np.diag(evals) @ evecs.T)[-1.83090639+0.j 2.87249557+0.j 4.0779703 +0.j]
------
[[ 0.78441158 0.52975091 0.32258712]
[ 0.47652876 -0.84764865 0.23326403]
[-0.39701237 0.02925297 0.91734696]]
True
Symmetric matrices
If A is symmetric you should use eigvalsh (H for Hermitian) instead: this is more robust and leverages the structures (you know they are real!)
Matrix operations
linalg.trace(A)# tracelinalg.det(A)# determinantlinalg.inv(A)# Inverse, consider NEVER using it though :)
Norms
print(linalg.norm(A, ord="fro")) # fro for Frobenius
print((np.sum(A ** 2)) ** 0.5)
print(linalg.norm(A, ord=2))
print((linalg.eigvalsh(A.T @ A) ** 0.5))
print(linalg.norm(A, ord=np.inf))5.313500804202211
5.313500804202211
4.0779703035175885
[1.83090639 2.87249557 4.0779703 ]
6.11483512481408
A = np.random.randn(3, 3)
print(linalg.norm(A, ord=np.inf))1.9635221329442685
Random generation, distributions, etc.
References:
- Good practices with numpy random number generators by Albert Thomas
- Numpy documentation on RandomState
- Random Widgets, by Joseph Salmon: Visualization of various popular distributions.
seed = 12345
rng = np.random.default_rng(seed) # can be called without a seed
rng.random()0.22733602246716966
Optimization
Goal: find functions minima or maxima
References:
- Scipy Lectures on mathematical optimization.
from scipy import optimizeFinding (local!) minima
def f(x):
return 4 * x ** 3 + (x - 2) ** 2 + x ** 4
def mf(x):
return -(4 * x ** 3 + (x - 2) ** 2 + x ** 4)
xs = np.linspace(-5, 3, 100)
plt.figure()
plt.plot(xs, f(xs))
plt.show()
Default solver for minimization/maximization: fmin_bfgs (see Wikipedia on BFGS)
x_min = optimize.fmin_bfgs(f, x0=-4)
x_max = optimize.fmin_bfgs(mf, x0=-2)
x_min2 = optimize.fmin_bfgs(f, x0=2)
plt.figure()
plt.plot(xs, f(xs))
plt.plot(x_min, f(x_min), "o", markersize=10, color="orange")
plt.plot(x_min2, f(x_min2), "o", markersize=10, color="red")
plt.plot(x_max, f(x_max), "|", markersize=20)
plt.show()Optimization terminated successfully.
Current function value: -3.506641
Iterations: 7
Function evaluations: 16
Gradient evaluations: 8
Optimization terminated successfully.
Current function value: -6.201654
Iterations: 5
Function evaluations: 12
Gradient evaluations: 6
Optimization terminated successfully.
Current function value: 2.804988
Iterations: 7
Function evaluations: 16
Gradient evaluations: 8

Find the zeros of a function
Find x such that f(x) = 0, with fsolve.
omega_c = 3.0
def f(omega):
return np.tan(2 * np.pi * omega) - omega_c / omega
x = np.linspace(1e-8, 3.2, 1000)
y = f(x)
# remove vertical lines when the function flips sign
mask = np.where(np.abs(y) > 50)
x[mask] = y[mask] = np.nan
plt.plot(x, y)
plt.plot([0, 3.3], [0, 0], "k")
plt.ylim(-5, 5)
optimize.fsolve(f, 0.72)
optimize.fsolve(f, 1.1)
optimize.fsolve(f, np.linspace(0.001, 3, 20))
np.unique(np.round(optimize.fsolve(f, np.linspace(0.2, 3, 20)), 3))
my_zeros = (
np.unique((optimize.fsolve(f, np.linspace(0.2, 3, 20)) * 1000).astype(int)) / 1000.0
)
plt.figure()
plt.plot(x, y, label="$f$")
plt.plot([0, 3.3], [0, 0], "k")
plt.plot(my_zeros, np.zeros(my_zeros.shape), "o", label="$x : f(x)=0$")
plt.legend()
plt.show()

Parameters estimation
from scipy.optimize import curve_fit
def f(x, a, b, c):
"""f(x) = a exp(-bx) + c."""
return a * np.exp(-b * x) + c
x = np.linspace(0, 4, 50)
y = f(x, 2.5, 1.3, 0.5) # true signal
yn = y + 0.2 * np.random.randn(len(x)) # noisy added
plt.figure()
plt.plot(x, yn, ".")
plt.plot(x, y, "k", label="$f$")
plt.legend()
plt.show()
(a, b, c), _ = curve_fit(f, x, yn)
print(a,"\n", b,"\n", c)
2.5945390998581153
1.1988808427139435
0.44817154460106
Displaying
plt.figure()
plt.plot(x, yn, ".", label="data")
plt.plot(x, y, "k", label="True $f$")
plt.plot(x, f(x, a, b, c), "--k", label="Estimated $\hat{f}$")
plt.legend()
plt.show()
Rem: for polynomial fitting, one can directly use numpy
x = np.linspace(0, 1, 10)
y = np.sin(x * np.pi / 2.0)
line = np.polyfit(x, y, deg=10)
plt.figure()
plt.plot(x, y, ".", label="data")
plt.plot(x, np.polyval(line, x), "k--", label="polynomial approximation")
plt.legend()
plt.show()/home/jsalmon/anaconda3/envs/peerannot/lib/python3.10/site-packages/IPython/core/interactiveshell.py:3505: RankWarning:
Polyfit may be poorly conditioned

Interpolation
from scipy.interpolate import interp1d, CubicSpline
def f(x):
return np.sin(x)
n = np.arange(0, 10)
x = np.linspace(0, 9, 100)
y_meas = f(n) + 0.1 * np.random.randn(len(n)) # add noise
y_real = f(x)
linear_interpolation = interp1d(n, y_meas)
y_interp1 = linear_interpolation(x)
cubic_interpolation = CubicSpline(n, y_meas)
y_interp2 = cubic_interpolation(x)
plt.figure()
plt.plot(n, y_meas, "bs", label="noisy data")
plt.plot(x, y_real, "k", lw=2, label="true function")
plt.plot(x, y_interp1, "r", label="linear interp")
plt.plot(x, y_interp2, "g", label="CubicSpline interp")
plt.legend(loc=3)
plt.show()
Images
from scipy import ndimage, datasets
img = datasets.face()
type(img), img.dtype, img.ndim, img.shape
print(2 ** 8) # uint8-> code sur 256 niveau.
n_1, n_2, n_3 = img.shape
np.unique(img)
# True image
plt.figure()
plt.imshow(img)
plt.axis("off")
plt.show()256

RGB decomposition
fig, ax = plt.subplots(3, 2)
fig.set_size_inches(7, 4.5)
n_1, n_2, n_3 = img.shape
ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256))
ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256))
ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256))
plt.tight_layout()
print(img.flags) # cannot edit...
img_compressed = img.copy()
img_compressed.setflags(write=1)
print(img_compressed.flags) # can edit now
img_compressed[img_compressed < 70] = 50
img_compressed[(img_compressed >= 70) & (img_compressed < 110)] = 100
img_compressed[(img_compressed >= 110) & (img_compressed < 180)] = 150
img_compressed[(img_compressed >= 180)] = 200
plt.figure()
plt.imshow(img_compressed, cmap=plt.cm.gray)
plt.axis("off")
plt.show() C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : False
WRITEABLE : False
ALIGNED : True
WRITEBACKIFCOPY : False
C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : True
WRITEABLE : True
ALIGNED : True
WRITEBACKIFCOPY : False

Convert a color image in grayscale
plt.figure()
plt.imshow(np.mean(img, axis=2), cmap=plt.cm.gray)
plt.show()
Blurring (fr: floutage)
img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.imshow(img_flou, cmap=plt.cm.gray)
ax.axis("off")
plt.show()Widget on blurring bandwidth (cannot be rendered online for the moment, only locally)
#| eval: false
%matplotlib widget
import ipywidgets as widgets
# set up plot
img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.axis("off")
plt.show()
@widgets.interact(sigma=(0.1, 200, 0.1))
def update(sigma=2):
"""Remove old lines from plot and plot new one"""
# [l.remove() for l in ax.lines]
img_flou = ndimage.gaussian_filter(img, sigma)
ax.imshow(img_flou, cmap=plt.cm.gray)Changing colors in an image
import pooch
import os
url = "https://upload.wikimedia.org/wikipedia/en/thumb/0/05/Flag_of_Brazil.svg/486px-Flag_of_Brazil.svg.png"
name_img =pooch.retrieve(url, known_hash=None)
img = (255 * plt.imread(name_img)).astype(int)
img = img.copy()
plt.figure()
plt.imshow(img[:, :, 2], cmap=plt.cm.gray)
fig, ax = plt.subplots(3, 2)
fig.set_size_inches(7, 4.5)
n_1, n_2, n_3 = img.shape
ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256), density=True)
ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256), density=True)
ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256), density=True)
plt.tight_layout()

References:
